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ABSTRACT 



The restoration of images focused on a photosensitive surface 
is treated from the standpoint of maximum-likelihood estimation, 
taking into account the Poisson distributions of the observed data, 
which are the numbers of photoelectrons from various elements of 
the surface. A detector of an image focused on such a surface 
utilizes a certain linear combination of those numbers as the 


optimum detection statistic. Methods for calculating the false- 


alarm and detection probabilities are proposed. It is shown that 
measuring noncommuting observables in an ideal quantum receiver cannot 


yield a lower Bayes cost than that attainable by a system measuring 


only commuting observables. I 
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Restoration of Images on Photosensitive Surfaces 


Methods for restoring images degraded by an optical system or by passage 
of the image-forming light through a turbulent medium often are based on a 
linear model of the optical process, in which the noise is independent of the 
desired image and combines with it additively. In these models it is necessary 
to postulate the statistical properties of the noise, which is usually assumed 
Gaussian, and to leave such parameters as its mean-square amplitude and its 
bandwidth to be measured separately. Such a treatment is incapable of deter- 
mining the fundamental limitations on restoring degraded images. 

Essaying to take into account the physical nature of the light and the 
statistical properties of the recording process, we have analyzed an imaging 
system in which the light from the object plane is focused on a photosensitive 
surface, from which it ejects photoelectrons. The surface is divided like a 
mosaic into a large number of small, insulated spots, from each of which the 
photocurrent can be measured. These measured values of the photocurrent con- 
stitute the data on which is to be based an estimate of the radiance of the 
object plane, which corresponds to the "true image". The numbers of photo- 
electrons emitted from the spots have Poisson distributions whose mean values 
are proportional to the instantaneous illuminance, which is the absolute square 
of the sum of the light fields created by object and background. These fields 
are spatio-temporal Gaussian random processes. The observed photocurrents 
are combined linearly to produce estimates of the object radiance, and the 
optimum linear combination determines the minimum mean-square error attainable 
by the restored image. 

In a paper, "Linear Restoration of Incoherently Radiating Objects", 
attached to this report, it is shown that the formula for the minimum mean- 
square error has the same form as that obtained in Wiener filtering theory. 
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but with the noise terms specified by the physical properties of the light 
fields and of the photosensitive surface. In particular, it is shown that 
under normal circumstances the shot noise arising from the stochastic nature 
of the photoelectron emissions far exceeds the noise associated with the 
random fluctuations of the object and background light fields. Image degrada- 
tion by diffraction at the aperture and by passage of the object light through 
a random phase screen, representing a turbulent medium, are both analyzed. 

As the shot noise predominates in causing random variations of the 
observed data, our knowledge of its properties can be made the basis of a 
technique for restoring images recorded by means of such a photosensitive 
mosaic. In developing the method, we describe the illuminance I(x) at point 
X of the image plane as the convolution of a true of geometrical image J(x) 
with an incoherent point-spread function (psf) , 


I(x) 


= J K(x - ^ 


y) J(y) d2y. 


( 1 ) 


The psf will be assumed normalized so that 
K(x) d^x = 1. 


/■ 


( 2 ) 


All integrals are carried out over the image plane, which will be taken finite 
for computational purposes. It is J(x) that we wish to estimate. 

The image plane is a photosensitive surface divided like a mosaic into a 
large number of small spots of area s, which emit photoelectrons when the light 
impinges on them. The product WT of the bandwidth W of the incident light and 
the observation time T is supposed to be so large, WT >> 1, that the numbers n^^ 
of photoelectrons from the spots are independent random variables with Poisson 
distributions. The mean value of n^^ is given by 


E(ni) = a s I(xi) , 


(3) 
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where is the center of the i-th spot and a is a constant proportional to 
the quantum efficiency of the surface and to the observation time and inversely 
proportional to the quantum hv of light energy. (More generally, a I(x) and 
a J(x) are averages over the spectrum of the incident light, weighted by 
q(v)/hv, where n(v) is the quantum efficiency at frequency v, the spectrum 
being assumed independent of x.) 

Let Iq = Jq be the uniform illuminance of the image plane in the absence 
of a scene; we then put 

J(x) = Jq + j (x), W 


suppose Jq known, and seek an estimate j (x) of the difference J(x) - Jq* The 
likelihood ratio of the numbers n^, i.e., the quotient of their joint probabil- 
ities in the presence and absence of an image, is then 


■ n[^] ' -p|- “ 

i * i 


[l(Xi) - Iq] 


!■ 


with 


I(x) - Iq = I K(x - y) j(;^) d^y. 

We seek estimates of j (x) at sampling points x^, and put j (x^^) 
I(xj^) = hi + Iq» approximating eq. (6) as 


(5) 


( 6 ) 

Si. 


Di = Cj. K^j = K(xi - Xj) s. (7) 

j 

We assume that the true image J(x) is a homogeneous spatial Gaussian random 
process with mean value Jg and covariance 


E[j(xp j(x2)] = 9(Xj - X2>, (8) 

and we form the matrix 9 whose elements are 

= E[5i Cj] = VCxi - Xj), (9) 
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and whose inverse Is u = |ivij^j||. The joint probability density function of the 
estimanda Ci is then 

' i j 

where M is a normalizing constant. 

The maximum-likelihood estimate of the true image is given by the set of 
values for which A{n^} is maximum. Equivalently, we maximize the 

logarithm 


) 


( 10 ) 


£n [A{nj^> pCfCi))] = = 


i i j 

where the are expressed in terms of the Cj by eq. (7). 
respect to 5 , we get 

- “ = •'«] ‘ “ • 
i m 


^ij ^i 

Differentiating with 


( 12 ) 


Solving these equations for we obtain 




m 


cp . K. . 

• Tm 1 T 




■Z.2L. "ij Llo + hi 
j 1 

■ Z - “ 


- a s I 


(13) 


where 



(14) 


If we define 

K'(x - y) 


J <p(x - u) K(u - y) d^u 


(15) 


and if we suppose the data nj> are given by nj^ = Zj^ s, where the z^ can be con- 
sidered as continuous random variables, = z(xj^), with conditional mean 
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(16) 


E(zi) = a Kg^), 


we can convert eq. (13) into an integral equation, 


j(x) 

= 1 

[*k'(x - y) 

z(y) 

lUy) ~ “J 

d2y. 

(17) 


J 

r 




i(y) 

- 1 

f K(y - w) j (w) d^w + 

lo- 

(18) 


J 


We can consider eq. (13) as a sampled version of eq. (17), with given in 
terms of Cj by eq. (7), which is the sampled version of eq. (18). It is then 
no longer necessary for the n^^'s to be integers. Equations (13) and (7) can 
be solved by iteration. 

The function F({5^}) in eq. (11) may have many minima, and there is a 
danger that the iteration method may settle into the wrong one. It is there- 
fore wise to start it as close to the absolute minimum as possible. One 
convenient starting point is the set of obtained from the linear least-squares 
estimate of the Cj's. Designating these by circumflexes, we put 


where 



(19) 


Ci = ~ Io)/as (20) 

j 

with the least-squares estimating filter obtained by solving the Wiener 
filtering equations* 


( 21 ) 

m m 

♦mj ■ <22) 

k,n 


* C. W. Helstrom, "Image Restoration by the Method of Least Squares", J. Opt. 
Soc. Am. vol. 57, 297-303 (March, 1967). 
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The last term in eq. (22) represents the noise, which is assigned a mean-square 
value equal to the variance of the numbers n^ when the illuminance is Iq* The 
least-squares estimate will be most nearly optimum when the contrast is low, so 
that I(x) = Iq, and when the numbers n^ of counts from each element of the 
surface are large enough, n^^ >> 1, so that the Poisson distribution can be 
approximated by a Gaussian with mean and variance given by eq. (3). Since 
the matrices K^j , , and have the Toeplitz form, finite Fourier transforms 

much simplify the computations involved in applying eqs. (19) - (22). 

In the iterative solution negative values of Jq + may be obtained. It 
is suggested that as these are unphysical, they be replaced by zero wherever 
they occur. There will then be no danger that I(y) or Iq + will vanish, 
since the kernel K(u) and the matrix elements kj^j are non-negative. If necessary, 
Jq + Cjjj could be constrained in the iteration to be at least equal to a positive 
number representing a uniform background illuminance. 
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False-Alarm and Detection Probabilities 
of Photoelectric Images 

1. False-Alarm Probability by the Method of Steepest Descent 


For the detection of an optical image providing an illuminance Ms(x) at 
point X of a photosensitive surface, in the presence of a background illuminance 
Mq(x), the optimum statistic has been shown to be* 


E Mj (3Ci)-| f 


(x) d^x, 


( 1 . 1 ) 


M 2 (x) = Mq(x) + Mg(x) , 


where x^ is the point from which the i-th photoelectron is emitted and t is a 
constant proportional to the quantum efficiency of the surface and such that 
the expected number n of photoelectrons during the observation interval is, in 
the presence of the image Mg(x), 



Integrals are taken over the entire photosensitive surface, and the sum in 
eq. (1.1) is taken over all the emitted photoelectrons. 

Here Hq will denote the hypothesis that the image is absent, background 
light alone being present; H 2 denotes the hypothesis that the image sought is 
present. The number n^ of photoelectrons from any small area 6A at point x^ 
of the surface has a Poisson distribution with mean values 


EtnilHo] 

= T Mg(Xi) 6A, 

(1.3) 

E[ni|Hi] 

= T M 2 (x^) 6A. 

(1.4) 


The statistic g is compared with a decision level gg; if g > gQ» the decision 

* C. W. Helstrom, "The Detection and Resolution of Optical Signals", Trans. 
IEEE, vol. IT-10, 275-287 (October, 1964). 
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for hypothesis Hi — image present — is made. More generally Mj (x) , j = 0, 1, is 
an average of the illuminance Jj(x, v) as a function of frequency v, weighted 
in accordance with the quantum efficiency n(v) of the photosensitive surface, 


Mj(x) 


j" n(v) Jj(x, v) dv/hv, j = 0, 1, 


where h is Planck's constant. Then t is the observation time. 
The false-alarm probability Qq is 


Qo = Pr{g > go|Ho], 


(1.5) 


and the detection probability is 
Qd = Pr{g > go |H i}. 


( 1 . 6 ) 


The probability distributions of g under the two hypotheses are generally 
difficult to calculate, and we shall endeavor to determine Qq and Qd from the 
Laplace transforms (moment-generating functions) 

fj(s) = E[e j = 0, 1, (1.7) 

of the distributions, which are given by 


fj(s) 


exp[uj(s)] 



f i 

Ml (x) 

expj 

T I 

Mq(x)_ 



0, 1; (1.8) 


these are related by 

f 1 (s) = f o(s - 1) . 


(1.9) 


The false-alarm and detection probabilities are then, by the inversion formula 
for Laplace transforms. 


Qo (so) 


^c+i“ 
J c — ±co 


^ J3/2.1 


( 1 . 10 ) 
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and 


•c+i° 


Qd(go) = 


1 - fi(sX „sgo 


e ds/2iTi. 


( 1 . 11 ) 


The contour of integration is a straight line parallel to and to the right of 
the imaginary axis in the complex s— plane. We concentrate here on the false- 
alarm probability; the technique for small values of the detection probability 
Qd is quite similar. 

Since the integrand of eq. (1.10) has no singularity at the origin or 
indeed anywhere in the complex plane, as an analysis of eq. (1.8) shows we 
can displace the contour to the left across the imaginary axis. Taking only 
the first term, 

/ 

s~^ ds/2Tri, 

we can complete the contour around the left half-plane and show that this 
integral vanishes. Hence the false-alarm probability is Qq = q(go)» where 


q(g) “ ■ I fo(s) e®® ds/2Tti, 

•'C 

the contour C lying parallel to and to the left of the imaginary axis. 


write this as 


q(g) “ I exp[po(s) + gs - Jin s] ds/2fTi 
'C 


I 


exp $(s) ds/2iri. 


where 


'i'(s) = Po(s) + gs - Zn(-s) 


( 1 . 12 ) 


We can 


(1.13) 


(1.14) 


is a complex phase. 

We apply the method of steepest descent by looking for the point s = sq 
at which the phase $(s) is stationary, that is, where 
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<I>'(s) = yo'(s)+g-s-^ = 0, S = so, 

denoting by primes differentiations with respect to s. The contour C is 
shifted until it passes vertically through Sq. Expanding the phase $(s) in 
a power series about sq, we write the integrand in eq. (1.13) as 


(1.15) 


exp $(s) = exp 


$(so) + Y (s - Sq)^ <I>"(so) 


+ (r!) ^ (s - sq)^ 'J’^^^Sq) 


r=3 


= exp 


[$(sq) + ^ (s - sq)^ <I>"(so)] 


OD 

1 (r!)"^ Cj.(s - Sq)^ 


(1.16) 


r=3 


where $^’^^(s) = d^ $(s)/ds’^, and where the coefficients c^. are derived from 
the series expansion of the exponential function. 


exp 


- 00 "1 _ _ 

(r!)-l (s - sq)’^ ‘t^’^^(so) = 1 +52 

-j -=3 - r=3 (1.17) 

Now we put for our variable of integration s = sq + iy, and integrating term 

by term on y we obtain 


q(g) = - [2 tt $"(sq)] ^ exp i3>(so) 


1 + 


ou 

E 


(- 1 ) 


'2k 


^ 2 ^ k! [$"(so)]’'j 


(1.18) 


In general, the series in brackets has only asymptotic validity and must be 
cut off when the terms, which at first decrease, begin to increase. The second 
derivative appearing here is, by eq. (1.15), 

<^"(so) = yo"(so) + so"^. 

A computer program to generate the coefficients Cj. in eq. (1.17) from 
these derivatives is simple to write. If in general we wish to take the 
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exponential function of a power series, 

^ (r!)“^ Cj. = exp ( ^ x®|, 

r=0 ' s=0 ' 

we can find the Cj-'s by the recurrence relation 

r-1 

Cr = ^ Vs ‘'s’ ‘'O = exp(bo), 

s=0 

which comes from Euler's formula for the (r-l)-th derivative of g'(x) exp g(x), 
where g(x) is the power series in the exponential function. Here bg = b^ = 
b 2 = 0, bj- = r > 3. 

fk') 

The derivative pg logarithm of the characteristic function 

(k) 

appearing in <I) (sg) can be written as 


f.YMl (s)\ 

k 

Ml (x) 



Mo(?)_ 


n-s 


pj^^s) = d^ pg(s)/ds^ = (-l)^TjMg(x) 

■ \--u '-W' / j L ' 

( 1 . 20 ) 

They will usually have to be integrated numerically, and a computer routine 
for generating all the necessary ones simultaneously will be expedient. 

The approximation of the false-alarm probability Qg derived here is 
closely related to the Chernoff bound, which is based on the inequality 


d'^x. 


exp si(g - gg) ^ U(g - gg). Si > 0, 


( 1 . 21 ) 


where U is the unit step function. Thus 

^00 

Qo = I U(g - go) Po(g) dg 


exp si(g - gg) pg(g) dg = e fg(-Si), 


( 1 . 22 ) 


* H. Chernoff, "A Measure of Asymptotic Efficiency for Tests of a Hypothesis 
Based on the Sum of Observations", Ann. Math. Stat. vol. 23, 493-507 
(December, 1952). 
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The value of is then picked to minimize the bounding function. Putting 
s = -sj , we must minimize 

£n fo(s) + sgQ = + SSOJ 

this requires 

Mo ' +80 = 0> (1.23) 

which except for the term -s“^ is the same as eq. (1.15). If sg' < 0 is the 
solution of eq. (1.23), the Chernoff bound asserts 


Qo(8o) - exp[po(so') + so'go]. (1-24) 

and the exponent here is nearly equal to 4>(sq'). In fact, since eq . (1.22) 
gives a bound for all positive values of s^, we can use instead of sg' the 
saddlepoint sg obtained from eq. (1.15) and assert the upper bound 


Qg(g) < exp[yg(sg) + Sgg], 


(1.25) 


which will be nearly as tight as that in eq. (1.24). The method of steepest 
descent thus yields both an approximation and a bound to the false-alarm 
probability . 

In determining the detectability of optical images, one fixes a false- 
alarm probability Qg and from it determines the decision level gg, which is 
then used in calculating the detection probability Q(j(8o). 1^ 1® therefore 
necessary to solve the pair of equations (1.15) and (1.18) for the value of g 
for which q(g) equals the pre-assigned false-alarm probability Qg . Here 
Newton's method is most expedient. The function q(g) can, through eq. (1.15), 
be treated as a function of Sg, which we designate by q(sg). Then if sg^ is a 
trial value of Sg, a new trial value is determined by the Newton equation 


so 


So 


+ 


Qo - q(so^) 
q'(sgl) ’ 


(1.26) 


where the prime denotes differentiation with respect to sg. In computing the 


I 
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derivative q'(so)> the series in the bracket in eq. (1.18) can be replaced by 
1, and the factor in front of the exponential function can be treated as 
constant; the rate of convergence to the solution will not be much affected. 
Thus we can use 

g'(so) = «>'(so) q(so), (1.27) 


where in the differentiation of iJ>(s) of eq. (1.14) g must be treated as a 
function of sq, 

4>'(so) = yo'(so) + g + g'so - so“^ = g'sQ (1.28) 

by eq. (1.15). On the other hand, eq. (1.15) gives 

Po"(so) + g' + = 0. (1.29) 

. Hence 

q’(so) = - q(so) [sQyo"(so) + Sq"^] (1.30) 


is to be used in the denominator of the Newton formula, eq. (1.26). 

In the beginning of the iterative calculation of sq, it will suffice to 
use for q(so) in eq. (1.26) the formula in eq. (1.18) with the bracketed series 
set equal to 1. At the end, as many terms of the series should be included as 
feasible. Rather than recalculating all the derivatives of ^>(s) at each new 
trial value of Sg through eq. (1.20), it may be more expedient to determine 
them from the Taylor series 


$(s) 


'y (r!) ^ ^^’^^Sg^) (S - Sg^)’^, 
r=0 


(1.31) 


where Sg^is a suitable trial value at some stage of the calculation. The 
coefficients in the bracketed series in eq. (1.18) will not need to be accurately 
evaluated if the contribution of the variable terms of the series is small. 
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2. The Probability of Detection 


For calculating the probability Q^j of detection, the method of steepest 
descent will be useful only for weak images, for which Qj is not much larger 
than Qo . The calculation will then be the same as what has just been described, 
except that fi(s) = fo(s - 1) is used in place of fo(s). Eventually, as 
increases, the series in the bracket in eq. (1.18) will cease to converge 
properly, and this method must be abandoned. 

For detection probabilities near 1, which are of most interest, an 
approximation method based on the analytical properties of fi(s) for large 
real values of s may be fruitful. In this analysis one may as well drop the 


constant last term 


Y = T ^ ^s(^) 

from the statistic g in eq. (1.1), using instead the modified statistic 


( 2 . 1 ) 


8l = g + Y, 


( 2 . 2 ) 


which is compared with the decision level g^Q = go + Y* The moment-generating 
function (m.g.f.) of g^ is 

fll(s) = E[exp(-sgi) | h^] = e fl(s). 

The probability density function (p.d.f.) of gj has a delta function at 
the origin, for gi = 0 when no photoelectrons at all are emitted during the 


observation interval. The probability of this event is 


Pl(0) 


r f 1 

= exp ~ T I (x) d^x 


(2.3) 


under hypothesis H^. The p.d.f. of gj, therefore, has the form 


Pl(gl) = Pl(0) <S(gi) + Pi'(gl) 
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The m.g.f. of the improper density function Pi ' (g) is 


/• 00 


fi(s) = 


PI ' (g) e dg = 


Pl(0) {exp 


fMo(x) 


s-1 


d^X 




(2.5) 


The behavior of Pi*(g) for small values of g, which are of principal concern 
when Qjj is near 1, is related to the behavior of f j^(s) for large, real values 
of s. The detection probability is given by 

rsio 

Pi '(g) dg 


Qd = 1 - Pi(0) - 


( 2 . 6 ) 


0 


and is therefore bounded above by 1 — Pj (0) independently of the decision level 
glO = go + Y* 

Let us consider, for example, a uniform background Mq(x) = Mq and a one- 


dimensional image of the form 

( A <p(x) , 0 < X < a, 

0 , X < 0, X > a, 


(2.7) 


with 9(0) = 0 and <p(x) a monotone increasing function of x. The image covers 
a range of width b in the y-direction. Since it is the behavior of Mg(x) in 
the regions where it is nearly zero that matters in what follows, this form 
is sufficiently general for the time being. 

The crucial term in eq. (2.5) is the integral in the exponent, which we 
write as 


J(s) = T / Mq(x) 
X = A/Mq . 


Mn (x) 


Ml (x) 


s-1 


d^x = T b Mr 


[1 + X cp(x)] dx. 


( 2 . 8 ) 


We wish to determine the dependence of this integral, and hence also of f^(x). 


on 


s when s >> 1. Let us take as an example cp(x) = x , v > 0. Then we get 
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J(s) = T b Mq 


ra 


0 

/•oo 


= T b M 


0 


(1 + X dx 


(1 + X x'^) dx; 


(2.9) 


the discrepancy introduced by extending the integral decreases exponentially 
with increasing s and can be neglected in this analysis. Evaluating the 
integral we get 


^OO 


-s.-l 


J(s) = T b Mq X 


V ^-1,, , .-(s-1) 

y (1 + y) dy 


= T b Mq v“i x“'' ^ r(v“^) r(s - 1 - v"^)/r(s - i) ( 2 . 10 ) 

in terms of the gamma-function F(x). If we now apply Stirling’s asymptotic 
formula for the gamma-function, we find for s >> 1 


J(s) ~ T b Mq v“^ X ^ r(v~^) s , s >> 1. 


( 2 . 11 ) 


Since J(s) << 1 when s >> 1, the m.g.f. fi(s) in eq. (2.5) is asymptotically 

fj^(s) ~ Pi(0) T b Mq v~^ X ^ r(v”^) s , s >> 1, (2.12) 

and its inverse Laplace transform is, for small values of gj , 


Pl'(gl) = Pl(0) T b Mq v~^ X 


-M~l v~^-l 

gl 


(2.13) 


The probability of detection is therefore approximately 

v-1 

Qd = 1 - Pi(0) - Pi(0) T b Mq (glQ/X)'^ . 


(2.14) 


We learn from this asymptotic analysis that the probability of detection 
depends strongly, when near 1, on the manner in which the image decreases to 
zero at its edge, which is reflected in the exponent v in eq. (2.9) and the 
exponent in eq. (2.14). If the image covers an infinite area, with Mg(x) 
decreasing to zero as |x| ^ a different asymptotic analysis is required. There 
will now in general be zero probability of zero counts, Pi(0) = 0, and the 


delta-function vanishes from the origin. 
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As an example we take a circular Gaussian image. 


Ms(x) ~ A exp(- X ^/a^). 


(2.15) 


to be detected against a uniform background Mq . The logarithm of the m.g.f. 
of gi becomes 

£n fii(s) = - T Mq [1 + a expC-r^/a^)] 



X |l - [1 + A exp(-r^/a^)] rdrdcp 


= - N I (1 + A e‘”) [1 - (1 + A e “) du, (2.16) 

'o 

where N = t Mg tt a^ is the mean number of counts under hypothesis Hg from a 
circle of radius a. In eq . (2.16) we make the substitution 


-u. 


V = Jln(l + A e ) 


(2.17) 


to obtain*, with b = £n(l + A), 


rb 


£n f 1 1 (s) = - N 


V/, -sv. 
e (1 - e ) 


1 - e 


-V 


dv 


= - N I (e + 1) dv 

'0 


r -V -(s 

- 5 ^ ^ 

Jo ^ 


-(s-l)v 


dv 


= - N (e - 1 + b) - N 


-V -(s-l)v 
e - e 


1 - e 


-V 


+ N 


-V -(s-l)v 
e - e 


1 - e 


-V 


dv 


= - N [A + £n(l + A)] - N [i|;(s - 1) + C] - N £n(l - e“^) - N AJ (s) 


= - N [A + £n A + i|j(s - 1) + C] - N AJ^(s), 


(2.18) 


where \p(z) is the logarithmic derivative of the gamma function, C = 0.577215 
is Euler's constant, and 


* I.S. Gradshteyn and I. W. Ryzhik, Tables of Integrals, Series, and Products 
(Academic Press, New York, 1965), §3.311, eq. (6), p. 304. 
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^ oo 


AJ^(s) = 


-(s-l)v 


1 - e 


-V 


dv. 


(2.19) 


Since over b<v<“, l-e^<l-e'^<l, we can bound AJ^(s) by 


-(s-l)b 


< AJ^(s) < 




-(s-l)b 


( 2 . 20 ) 


s - 1 ■' \1 - e -/ (s - 1) 

and we see that it decreases to zero exponentially as s Hence asymptotically 

fll(s) ~ exp{- N [X + iln X + C + ijj(s - 1)]}, s >> 1. (2.21) 

By using the asymptotic expansion* of the function 4 ^( 2 ), 

^ no , 1 (2.22) 


'li(z) 


£n z - 1- 0(z ^) , R£ z >> 1, 

2 . 7 . 


we obtain 


flj(s) ~ exp{- N [X + C + £n(Xs) - + 0(s~^)]} 


2s 

= (X.r® U-||+0(s-2)]. 

The p.d.f. of the statistic is therefore, for gj << N, 

. . -N -N(X+C) N-1 3 . r./- 2 m 

Pl(gl) ~ ^ e r(N) ^ gj [1 - j gl + 0(gi"^)], 

and the probability of detection is 


,-N -N(X+C) . ,^,-1 N ,, 3N . 2s s 

Qd ~ 1 - X e [r(N + 1)] ^ gio (1 - 2) 0 ( 810 ) ) 

where giois the decision level on the statistic gj . 

For detection probabilities in the intermediate range where neither of 
the asymptotic methods just described will work, one can calculate the characteristic 
function E(e^*^®l |Hi) = fi(-io)) of the distribution and take its inverse Fourier 
transform 

Pl(Sl) = 


fl(-iu) e 


-iu)g] 


du/2ir 


* A. Erddlyi > Higher Mathematical Functions (McGraw-Hill, New York, 

1953), vol. 1, §1.18, eq. (7), p. 47. 
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by one of the fast computer algorithms now available. This can then be inte- 
grated numerically to give the detection probability. Alternatively, the Gram- 

Charlier series can be tried, although it appears to be successful only for 

A 

unimodal distributions that closely resemble the Gaussian. 


C. W. Helstrom, Statistical Theory of Signal Detection (Pergamon Press, 


Oxford, 2nd ed., 1968), pp. 219-222. 
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Quantum Detection Theory 


An eighty-page review article on quantum detection theory has been 
prepared for publication in Progress in Optics , E. Wolf, editor. 

Quantum detection is a form of statistical hypothesis testing adapted 
to quantum-mechanical rather than classical laws of statistics. It has been 
formulated in terms of the conventional rules of quantum mechanics, one of 
which is that only observables associated with commuting operators can be 
measured simultaneously on the same system. The possibility of simultaneously 
measuring noncommuting observables has been studied in recent years, and it 
has been suggested that this would permit quantum receivers to attain lower 
average error probabilities or Bayes costs. In a paper attached to this 
report it is shown that no reduction in the minimum Bayes costs can be 
achieved in this manner. Measuring noncommuting observables requires an auxiliary 
apparatus that initially possesses no information about the state of the 
quantum receiver, and there is no way in which introducing such an apparatus 
can lead to a lower Bayes cost than the minimum attainable by measuring 
commuting observables on the receiver alone. Quantum-mechanical Cramfir-Rao 
bounds on the mean square errors of unbiased estimates of parameters of a 
signal in a quantum receiver can likewise not be lowered by measuring noncom- 
muting observables. 
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Noncommuting Observables in 
Quantum Detection and Estimation Theory 

Carl W. Helstrom* 

Department of Applied Physics and Information Science 
University of California, San Diego 
La Jolla, California 92037 


Abstract 

In quantum detection theory the optimum detection operators 
must commute; admitting simultaneous approximate measurement of 
noncommuting observables cannot yield a lower Bayes cost. The lower 
bounds on mean square errors of parameter estimates predicted by 
the quantum-mechanical Cramer-Rao inequality can also not be 
reduced by such means. 


* This research was supported by Grant NGL 05-009-079 from the 
National Aeronautics and Space Administration. 
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Quantum detection and estimation theory has been developed within the 
conventional framework of quantum mechanics, one of the principal tenets of 
which is that only observables associated with commuting operators can be 

[1-3] 

simultaneously measured on the same system. It has been suggested that 

this formulation is too restrictive, that noncoramuting operators can be at 
least approximately measured on the same system, and that to include this 
possibility may permit more effective detection, as measured by a lower average 
Bayes cost.^^’^^ We wish to show that no such improvement can be expected. 

The simultaneous measurement of noncommuting observables has been 
treated by Gordon and Louisell.^^^ In order to approximately measure certain 
such observables on a quani um-mechanical system S, it is made to interact for 
a time with a second system A, termed the apparatus. It was shown that a 
suitably defined ideal measurement yielding approximate values of the noncom- 
muting observables can be based on the outcome of measurements of commuting 
observables on the apparatus A, or more generally on both S and A. What we 
must therefore do is apply quantum detection theory — with its restriction to 
commuting observables — to the combined system S + A. 

Suppose we are to decide among M hypotheses Under 

hypothesis Hj the density operator for the combined system at time t is 
S+A 

Pj (t) in the Schrddinger picture. If at an earlier time tg the density oper- 
. S+A /S-L 

ator IS pj (to)> the two operators are related by 

p^'^^(t) = U(t, tg) tg), (1) 


with 


rt 


U(t, tg) = exp 


-1 


Hdt'/fi 


( 2 ) 


Jtg 


where H is the Hamiltonian operator for the combined system S+A and fi is 
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Planck's constant h/2-rr. The operator U is unitary; that is, with U its 
Hermitian adjoint, UU^ equals the identity operator 1. 

Let {Hj} be a set of commuting projection operators forming an M-fold 
resolution of the identity. 


M 



(3) 


j=l 

On the combination S + A we are to measure these M projection operators at time 
t, and if the k-th yields the value 1, hypothesis is selected as true.^^^ 

The average cost is then 


C 



Tr[p^^+^(t) n.]. 


( 4 ) 


i=l j=l 

where is the prior probability of hypothesis Hj and is the cost of 
choosing when Hj is true. Let {llj(t)} be the projection operators that 
minimize C when the system S + A is observed at time t; we call these optimum. 
Then by (1) the operators 


Hj(tQ) = U”^(t, tg) HjCt) U(t, tg) (5) 

f 

will minimize C when S + A is observed at time tg. Because of the unitarity 
of U(t, tg), the set {llj(tg)} also forms an M-fold resolution of the identity 
into commuting projection operators, and the Ilj(tg) are optimum at time tg. 
Since the minimization is carried out over all possible M-fold resolutions of 
identity, the minimum Bayes cost must be independent of the observation 

time t. 

Now let us roll time back to an epoch tg before the system S has come 
into contact with the apparatus A. In the Schrbdinger picture this amounts to 
applying the inverse unitary transformation u"*'(t, tg) to the state vectors of 
the combined system S + A. Because S and A are independent at this time to. 
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the density operators ^ must now have the factored form 
pj^^(to) = Pj (to) P (^o)» j “ 


( 6 ) 


Furthermore, as the apparatus A before the interaction has no information 
about which hypothesis is true, p^(to) in (6) must be independent of j. 
The Bayes cost is now 


C 





P.^to) n. 


(7) 


L i=l j = l -J 

Since S and A are completely uncoupled, and the state of A is independent 

of which hypothesis H. is true, there is nothing to be gained by observing A. 

S A A . 

The optimum projection operators Ilj^(to) factor as 11^^ (to) 1 > where 1 is the 

S 

identity operator for the apparatus A, and the set {ITj (to)) forms an M-fold 

s 

resolution of the identity 1 for the system S, minimizing the Bayes cost 


M M 


C 


S 




Tr[pj^to) 




( 8 ) 


i=l j=l 

Since Tr[p'^ 1^] = 1, the minimum value of is also the minimum value of C 
in (7) and equals the time-independent minimum Bayes cost The decision 

among the M hypotheses made at time to is based entirely on the measurement 
of commuting observables on system S. 

Similar considerations apply to estimating the m parameters 0 = 

O 

(0^, 02 ,..., 0ju) of the density operator p (0) of a quantum-mechanical system 

S. A version of the Cramer-Rao inequality sets lower bounds to mean square 

[3] 

errors of unbiased estimates of 0i, 02,--.> 6,jj» Let be an operator 

A 

whose measurement on S yields an unbiased estimate 0j of the j-th parameter; 

0j must be an eigenvalue of . Although in order to be measured simultaneously 
on the same system the operators Xj must commute, the analysis leading to the 
lower bounds given in [3] does not require commutativity of the operators X j . 
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For the class of aomrmting operators yielding unbiased estimates of the para- 
meters 0 there will exist lower bounds on the mean square errors, and those will 
be greater than or at least equal to the bounds derived in [3] . 

Again including the possibility of measuring noncommuting observables 
cannot lead to lower bounds smaller than those in [3] . In order to measure 
such operators even approximately, a measuring apparatus A must be allowed to 
interact with the system S, and according to Gordon and Louisell's treatment 

of the process, commuting operators will at the end be measured on the combined 

r 6 1 S+A 

system S + A. In the Schrbdinger picture the density operator p (6, t) 

for S+A will have a time dependence similar to that in (1) . 

Referring to (7) of [3] we see that the symmetrized logarithmic deriva- 
tives (SLD) Lj (t) appropriate for determining the Cramer-Rao lower bounds when 
the measurements of Xj are made at time t are related to those appropriate for 
measurements made at tg by 

Lj(t) = U(t, to) Lj(to) U'*’(t, to). (9) 

Then (13) of [3] shows that the matrix A that sets the lower bounds is indepen- 
dent of the time t of observation, again because of the unitarity of the 
operator U(t , tg) • 

Once more we move back to an epoch tg before the system and the apparatus 

have interacted. The density operator p^'*'^(0, tg) factors as p^(6, tg) p^(tg), 

A 

where the density operator p (tg) of the apparatus A is independent of the 
estimanda 0. The SLD operators for calculating the lower bounds are now the 
solutions of the operator equations 

8p^(0, tg)/30^ = I [p^(0, tg) P^(0, tg)], (10) 

and they act only on system S, commuting with p and all other operators on the 
apparatus A. When taking the trace over the states of A to form the elements 
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of the matrix A, the density operator is replaced by 1, and the lower 
bounds depend only on p (0, tg)- Thus the apparatus A cannot help estimate the 
parameters 6 of S with smaller mean square errors than the lower bounds cal- 
culated by the quantum-mechanical Cramer-Rao inequality as applied to the 
density operator of system S alone. 

In [3, p. 238] lower bounds were calculated for unbiased estimates m^^ and 
niy of the components of the complex amplitude y = nij^ + imy of a simple harmonic 
oscillator, which might represent a mode of the field in an ideal receiver in 
the presence of thermal noise. Those bounds are 

Var (N + i) » Var Ay > (N + y) . 

where N is the mean number of noise photons. The noncommutativity of the SLD's 
Ljj and Ly used to derive these bounds does not invalidate them. It can be 
shown that if the mode is coupled with an ideal amplifier whose gain is high 
enough to raise the oscillator variables to the classical domain where they 
commute, error variances 

Var m^ = Var Ay = y (N + 1) 

can be attained [8]. It is unknown whether commuting operators can be found 
whose measurement will yield unbiased estimates Aj^ and Ay with variances lying 
between y (N + y) and ^ (N + 1) • 

The measurements we need to make on a quantum-mechanical system S for 
testing hypotheses or estimating parameters will always have to be effected by 
means of an auxiliary apparatus A, and this apparatus, subject to thermal and 
quantum fluctuations of its own, will ordinarily introduce additional random 
uncertainties. Each measurement procedure will have to be analyzed to determine 
what error costs it entails. Detection theory and estimation theory seek 
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lower bounds on these costs, and in doing so they minimize with respect to the 
entire class of possible detection or estimation operators that can be applied 
to the system. The resulting bounds are independent of the time of observation, 
and they cannot be reduced by using any auxiliary apparatus that initially 
posseses no information about the state of the system. 
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Linear Restoration of Incoherently 


Radiating Objects 


Carl W. Helstrom* 
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ABSTRACT 


Light from an incoherently radiating object and background 
light are focused onto a photosensitive mosaic, the currents from 
whose elements constitute the data on which is based a least- 
squares linear estimate of the radiance at points in the object. 

By comparison of the mean square error with that given by Wiener 
filtering theory, the equivalent noise spectral density for use in 
the latter is shown to consist of a shot-noise term and a term due 
to the random fluctuations of the incoherent light; the former 
predominates under most circumstances. Turbulent distortion of 
the image after passage of the rays through a random phase screen 
is also treated from this standpoint. 

* This research was supported by Grant NGL 05-009-079 from the 
National Aeronautics and Space Administration. 
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The restoration of degraded images is often based on a linear model of 
image formation: a two-dimensional object function is pictured as having 

passed through a spatial filter, to whose output spatial random noise has 
been added; the sum of these constitutes the observed data.^ ^ In the 
approach known as Wiener filtering, the data are passed through a linear 
spatial filter whose transfer function is designed to reproduce the original 
object function within the smallest possible mean square error. In 
actuality, most objects radiate incoherently, and it is not the light field 
they create that is wanted; that field is a complex spatio-temporal gaussian 
random process, whose instantaneous values are of no interest. Rather the 
observer wants to know the distribution of radiance across the object, and 
the radiance is related to the mutual coherence function, which is the statis- 
tical average of a quadratic functional of the light field. Furthermore, the 
linear filtering theory does not explicitly distinguish the two basic types 
of noise, that due to the inherent fluctuations of the light fields of object 
and background, and that associated with the process of recording the light. 
Instead it makes certain assumptions about the spatial spectrum of the noise 
and requires its strength to be measured separately. This theory is inadequate 
to determine the basic limitations set by nature on the restorability of 
degraded images. 

A previous paper showed how the statistical properties of the information- 
bearing quadratic functional of the light field modify the usual Shannon 

Q 

formula for the information transfer from an incoherently radiating object. 

Here we shall treat the linear processing of the image illuminance that seeks 
an estimate of the radiance distribution in a distant object plane. In order 
to incorporate the random behavior of the recording medium, we shall assume 
that the light from the object is focused on a photosensitive surface divided 
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like a mosaic into many small, contiguous regions, the photoelectronic currents 
from each of which make up the observed data. These data are combined 
linearly to provide estimates of the radiance at discrete points of the object 
plane. 

Only spatial filtering of this kind will be considered. Temporal 
filtering for reducing the estimation error on the basis of color differences 
between object and background will be disregarded. The object light and the 
background light will be assumed to pass through a narrowband filter, over 
whose passband both have uniform spectral densities. 

We begin by working out the means and covariances of the photocurrents 
in terms of the mutual coherence function of the incident light, which in 

I 

turn is related to the radiance distribution of the object plane, our estimandum, 
and to the background radiance. The object radiance is considered as one of 
an ensemble of spatial random processes, which form the class of objects to be 
examined. The squared error of the radiance estimate is averaged over this 
ensemble of object processes, as well as with respect to the statistical dis- 
tributions of the light fields and the photoelectrons. The spatial filtering 
that minimizes the resultant mean square error is then easily determined. 

The final expressions for the mean square error will be found to have 
the same form as those derived from the Wiener filtering theory, and by 
comparing them we can determine the spatial spectral density that must be 
assigned to the random noise in utilizing the Wiener technique. This equi- 
valent noise density is the sum of a constant term due to the randomness of 
the photoelectron emissions ("shot noise") and a spatially variable term 
arising from the fluctuations of the light fields produced by object and 
background. Finally, since removing distortion by atmospheric turbulence is 
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a primary goal of image restoration, the effect of a simple kind of turbulence 
the random phase screen — will also be analyzed in this way. 
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1. The Photocurrents 

The sysL \n to be analyzed is shown in Fig. 1, The object plane 0 emits 
incoherent light, which is received at the aperture A of an observing optical 
instrument, A narrowband filter, not shown, cuts out background light outside 
the temporal spectrum of the object light. The aperture A contains a lens 
to focus the light on the image plane I, which consists of a mosaic of isolated, 
but close photosensitive spots of area s. The light ejects photoelectrons from 
these spots, and the resulting integrated currents are measured and consti- 
tute the data on which estimates of the object radiance are to be based. The 
spots are so small and so close together that ultimately we shall treat them 
as infinitesimal. 

The object pJ.ane is likewise sampled at a dense set of points uj^^ = 

(^kx* radiance estimated as a linear combination of 

the integrated currents Ijjj, 

B(yk) ^km^m const., (1.1) 

m 

where the weights Lj^ introduced by the processing matrix L are selected to 
minimize a mean square error; the constant has a known value. 

Dividing time into small intervals At, we write the current from the m-th 
spot as 

Im(t) (1-2) 

i 

where n^^^Ct^) is the number of photoelectrons emitted during the i-th interval, 
and f(t) is the current pulse resulting from each emission; At is much smaller 
than the duration of f(t). For a given light field at the image plane, the 
numbers are random variables with a Poisson distribution, and they are 
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statistically independent from one spot to another and from one interval to 
another . 

The object, aperture, and image planes are assumed to be so far apart, 
and the aperture so small, that all light rays can be considered paraxial. If 
for simplicity we suppose the light linearly polarized, it suffices to use a 

scalar theory, and we represent the field at point x of the image plane by the 

analytic signal t|;(x, t) e where = 2 t ^ c /\ is its central angular frequency, 

A is its wave'. ngth, and <l^(x, t) is the complex envelope. The conditional 

expected value of the number is given in terms of this field by 

E[nm(ti)|4'] = -| t^)|2 At, (1.3) 

where a = n/flfi is the quantum efficiency n of the surface divided by the 
quantum hfi of energy of the light. The conditional mean of the current Ijjj 
from the spot centered at Xjjj is then 



The current from each spot will be integrated by a filter with a long 
time constant T in order to determine the total charge ejected by the light 
during the observation interval. This filter, acting on each component pulse 
f(t) of the current, replaces it in Eq. (1.2) by a new pulse whose duration is 


of the order of the integration time T. We can therefore regard the data as 
represented by I^jj(t) of Eq. (1.2), but with a pulse shape f(t) that is positive 
and of duration T. If we suppose f(0) = f(t), a useful definition of the 

duration T is 
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[f(t)]2 dt/[f(0)]2; 


(1.5) 


T is in effect the observation time. Without loss of generality we put f(0) = 1. 

In order to work out the mean square error of our estimate, we shall 

need the covariance of the currents from two different spots. It can be evalu- 

9 10 

ated as in the semi-classical analysis of the Hanbury-Brown-Twiss effect, ’ 
which utilizes the equality of the variance of the Poisson distribution of 
nm(tj^) with the mean value given in Eq. (1.4), and which also involves taking 
the expected value of a product of four correlated complex gaussian variables 
by the formula^^ 

The resulting covariance is 


{Iro(t), Ip(t)} = <Im(t) Ip(t)> - <Im(t))<Ip(t)> = 


l'Fi(Xjjj, Xp)|2 


|x(ti - T2) p f(t - Ti) f(t - T 2 ) dTidx2 


+ i'i(Xj^, ^) 6 


[f(t - t)] 2 dx. 


( 1 . 6 ) 


where 6jjj is the Kronecker delta and Tj(xi, X 2 ) is the spatial part of the 


mutual coherence function of the light field at the image plane, 
■|<'P(5l. t^) t^)} = 'l'i(?i. 52^ “ ^ 2 )' 


(1.7) 


The temporal part x(tl “ ^ 2 ) is the Fourier transform of the absolute square 
X(uj) of the transfer function of the input filter and is so normalized that 


X(0) = 1, 


x(x) = I X(m) e dio/Zn. 


Here angular frequency m is defined with respect to Q as origin. The angular 
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brackets in Eq. (1.6) refer to averages with respect to the distributions of 
both the light fields and the photoelectron emissions. 

If we define the bandwidth W of the filtered light as 


r 

/*co 

1 |x(t) |2 dx = 

1 X(m) dm/2Tr 

/-OO 

^ —00 


W = |x(0)|y I |x(t)| 2 dx = I I X(w) dm/2Tr| j j |x(m)|2 dm/27r, 

( 1 . 8 ) 

and if we assume WT >> 1, as is normal in practice, the double integral in 


Eq . (1.6) can be simplified, and the covariance of the currents becomes 


{Im(t), Ip(t)} = a2s2 [^^(x^, Xp)|2 T/W + as ^) 6^^ T.(1.9) 

In terms of the mutual coherence function, the mean value of the current from 
the m-th spot is, by Eq. (1.4), 

<Im(t)> = as Vj(x^, x^) pT, (1.10) 

where p is a positive constant defined by 


P 


r-1 


f(t) dt/f(0) 


( 1 . 11 ) 


For rectangular current pulses p = 1. 

We must now evaluate the spatial coherence function l'i(x^, x^) in terms 
of the radiance distribution B(u) of the object plane and the spectral density 
N of the background light. 
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2. The Mutual Coherence Function 


The field t) at point x of the image plane is given in terms of 

the aperture field t) by 

where S^Cx, r) is the point-spread function (psf) between aperture and image 
plane. In the paraxial approximation 

^ 1 ^ 5 ’ = i(ARi)“^ exp 

where k = n/c = 2 tt/A is the propagation constant and is the distance from 
the aperture to the image plane. The last term in the exponent is the phase 
shift introduced by the lens, whose focal length is F and which focuses the 
object plane at distance R onto the image plane, 

F-1 = R-1 + R^-1. (2.2) 


IkRi 4- f- 


r - X 


ik o 
— r 
2F ~ 




( 2 . 1 ) 


The aperture field consists of two statistically independent parts, the 
light from the object and the background light, 

'('a(^» ^ f y) 'f'o^y» + 'l'n^5» (2.3) 

Jo 

where \po(u, t) is the field at point u of the object plane, 82 ( 5 , u) is the 
psf between object and aperture, and '/'„(£> t) is the background component of 
the aperture field. In the paraxial approximation 

ik 


S 2 (r, u) = i(AR)~^ exp ikR + |u - rj^ . 


(2.4) 


As the background light is spatially incoherent, its mutual coherence 
function after filtering has the form 




(2.5) 


9 



where 6(r) is the two-dimensional delta-function. The object light is spatially 
incoherent at the object plane, where its mutual coherence function is 

■|<)|Jq(u^, tj) iJJq*(u 2, t^)) = (x2/4tt) B(u^) (S(u^ - U 2 ) x(ti - t2> 

( 2 . 6 ) 

if we suppose the temporal filtering already applied. Here B(u) is the radiance 
of the object plane integrated over the passband of our input filter. It is 
decomposed as 

B(u) = Bq + b(u) . (2.7) 

into a known average level Bq and a deviation b(u). It is b(u) that the 
system is to estimate; it represents the scene of interest. 

When we combine all these equations, we find for the spatial coherence 
function of the light at the image plane 


X 2 ) = (AttRi^) 1 exp [ik(x^2 _ 


X <(Bo' 


exp [ikf(x 2 - x^)/Rj] d^r + 



where 


(AR)~2 / / / b(u) exp 

'A A *2o 

Bq' = Bq + 4^N/a2 


ik , . , ik , . 

— (X2*r2 - xj-rj) + — u- - r^ ) 


X d2ud2rjd2r2 /> 


( 2 . 8 ) 


(2.9) 


and A and 0 indicate integrations over the aperture and the object planes, 
respectively. The uniform radiance level Bq adds to the background density, 
and we merge the two into a single uniform radiance level Bq'. In our discussions 
4 ttN/A^ will be neglected in comparison with Bq; it can easily be restored by 
using Eq. (2.9) . 

It is convenient to rescale the coordinates in the image plane so that 
the geometrical image of a point u in the object plane also has coordinates u; 
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we put 


V = - R x/R^ 


( 2 . 10 ) 


and write the spatial coherence function as 


~ 2 ^ A( 47 tRi^) ^ exp [ikRjCvj^ - V2^)/2 Rq^] 


Bo'.y(v^ - Y2) + (AR) ^ A I b(u) - u).j^*(v2 - u) d^uL 


= A 


= A-1 


I^(r) exp (iku*r/R) d^r 


is the Fourier transform of the indicator function of the aperture, defined 




( 2 . 13 ) 
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3. The Radiance Estimate 


If we now use Eqs. (1.10) and (2.11) to evaluate the expected value of 
the curreiit from the m-th spot on the image plane, we obtain 


<Im> = If 


+ p s 


b(u) |.i^(Vjjj - u) p d^u, 


(3.1) 


Iq = ctpTs A Bq' (4tiRi2)-1^ 

Cl = aT A2(XR)-2(ATrRj2)-l^ (3 2) 


where v^ is the scaled coordinate vector of the center of the m-th spot. 

Except for the known constant Iq, the mean value of Ijjj is proportional to the 
weighted integral of the true radiance deviation b(u) over the neighborhood of 
the object point Vj^. The kernel \^i.v - u) j 2 is the incoherent point-spread 
function (psf); its diameter is of the order of the conventional resolution 
interval 6 = XR/a, where a is the radius of the aperture, taken as circular. 

If b(u) varies only slightly over a distance 6, 1,^ - Iq itself provides 
a good estimate of b(vjjj); the linear processing of the currents 1^^^ is intended 
to improve that estimate, and the new estimate of the radiance deviation at 
point Ujjj of the object plane is taken as 



(3.3) 


n 

with the coefficients L to be determined as shown hereafter. 

mn 

The expected value of the radiance estimate b(Vjjj) can be obtained from 
Eqs. (3.1) and (3.3), 


= P CiS / b(y) |d?(Ujjj - v) 

n *^0 

where s is the area of an element of the image plane. 


^ d2y, (3.4) 

We now assume that s is 
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so small, and the spots so close together, that this summation can be approxi- 


mated by an integral. 


<b(u)> = p Ci(Ri/R) 2 / / L(u, w) l^(w - v)|2 b(v) d^wd^v, (3.5) 


where we introduce a smooth weighting function L(u, w) whose values at the 
sampling points are L(uu,, u^^) = The factor (R^/R)^ arises from the 

scaling adopted in Eq. (2.10). 
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4. The Mean Square Error 


Any system for image restoration must be evaluated with respect to a 
definite class of objects. The statistical properties assigned to the class 
reflect the nature of the objects whose images the system is designed to 
restore. Here the objects are regarded as having radiance distributions 
B(u) = Bq + b(u) that are homogeneous two-dimensional random processes of 
mean value 

E B(u) = Bq (^-1) 

and autocovariance function 

E[b(u) b(v)] = cp(u - y), Og2 = cp(0). (4.2) 

The width of cp(u) as a function of |u| represents the size of typical details 
in the object plane; the ratio Og^/B'Q^ specifies the mean square contrast. 

The Fourier transform of cp(u) provides the distribution of spatial frequencies 
in objects of the class. Objects in which sharp edges predominate, for 
example, will have spatial spectra decreasing slowly at high spatial frequencies. 

As the actual object radiance is unknown a pTiori, the system can be 
designed only to attain a certain average performance with respect to the given 
class of objects. If it does so, it can be expected to restore efficiently 
most objects of the class. The measure of performance adopted here is the 
common mean square error 

S = E([b(u) - b(u)]^), (A. 3) 

in which the angular brackets denote, as before, averages over the distributions 
of the fields and the photoelectron emissions, and the symbol E now denotes an 
expected value over the class of anticipated objects. The set of coefficients 
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Ljnn> the weighting function L(u, v) of which they are samples, is to be 
chosen to minimize 

The mean square error can be written as 

^ = E<[b(u) - <b(u)> + <b(u)> - b(u)]2> 

= E Var b(u) + E[<b(u)> - b(u)]2, (4. A) 

where Var stands for the variance with respect to the field and photoelectron 
distributions. From Eqs. (3.3) and (1.9) 

MSm) ll»p (I„. Ip> = 

n p 


{«^s2|4<j(x„, Xp)|2 T/W + 

n p 

^mn 5n> ^ 


n 


a2(T/W) (Ri/R)** / / L(^, L(um» Y 2 ^ Y 2 ^l^ ^^Yi^^^Y 2 


aT(Rj/R)2 / [L(ujn, v)]2 V^(v, v) d2y 


(4.5) 


in the limit s ->■ 0, where Y 2 ^ spatial coherence function given 

by Eq. (2.11). Substituting from Eqs. (3.5) and (4.5) into Eq. (4.3), using 
Eq. (2.11), and averaging by means of Eq. (4.2) over the ensemble of radiance 
deviations b(u), we obtain for the mean square error 


a2(XR) 


u, Vj) L(u, y^) |Bo^|.i^(y^ ~ Y 2 ^l^ 

- Uj) J?(uj - y2) ^(Y2 ~ V2^ '^^^2 ~ Yi) 



X d2ujd2y2 > d2vjd2y2 


+ C2 Bo'(AR)2 A-1 / |L(u, v)]2 d2y + 
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P^C 2 ^ I L(u, v^) L(u, Y 2 ) - Wi)|2 I^(Y 2 “ ^2^1^ “ ~2^ 


X d^v^d^V2d2wjd^W2 


- 2 p C 2 / L(u, v) |^(v - w) ^ cp(w - u) d^vd^w + cp(0), 


(4.6) 


with C 2 = Cl (R i/R)2. 

Because of the homogeneity of the object and of the optical imaging in 
our paraxial approximation, the weighting function L(u, y) that minimizes S 
will depend on u and y only through u - y. We introduce its spatial Fourier 
transform A(r), defined by 


L(u, v) = (AR) £ 2 ^ I A(r) exp [ikr»(u - y)/R] d^r; (4.7) 


the transform variable r is scaled so that it matches the coordinates in the 
aperture plane. Similarly we define the spatial spectral density of the 
radiance deviation b(u) through the Fourier integral 


cp(u) = / o(r) exp (ikyr/R) d^r. 


(4.8) 


In terms of these the mean square error can be written as 


S’ = (WT) 


/ |A(r)|2 B 52 I^^^(r) + / s) $(s) d^s d^r/A + 


C I |A(r)|2 d^r/A + / |pA(r) I, ' (r) - Ip <i>(r) d^r. 


(4.9) 


where C = 4tr Bo'/A^aT, 


ir (r) = / Ia(§) Ia(§ - V d^s/A 


(4.10) 


is the self-convolution of the indicator function of the aperture, and 


^A'~^ -A' 


is a quadruple convolution. For a circular aperture of radius a, both (r) 
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and s) vanish outside a circle of radius 2a. Their maximum values occur at 

r = s = 0 and equal 1. The region where 0 is called the convolved aperture. 

The expression in Eq. (4.9) has the same form as the mean square error in 

( 2 ) 

Wiener filtering; 4>(r) is the object spectral density, p (r) the effective 

optical transfer function, and 


$^(r) = 4 tt Bo'/A^aTA + 


A“1 (WT) 


-1 


+ I s) Ks) d^s 


(4) 


(4.12) 


the equivalent spatial spectral density of the noise. In 3>j^(r) the first term 
represents the shot noise due to the photoelectrons, and the second arises from 
the fluctuations of the light field. Thus the mean square error is 

J lA(r)|^ '^n^5^ j 1pA(h) ^ ~ (^*13) 

in which the first integral gives the contribution of the noise and the second 
represents the mean square bias of the estimate, averaged over the ensemble of 
radiance distributions. 

The filter transfer function A(r) that minimizes the mean square error 
is determined as previously described. As usual in linear processing, spatial 
frequencies in the object that when multiplied by XR fall outside the convolved aper- 
ture are not restored to the image. When the aperture is so large that diffraction 
can be neglected, only the shot noise and the fluctuations of the light field 
prevent perfect image restoration. The equivalent noise density is then constant, 
and the minimum mean square error is 


S. 


[1 + H $(r)/$(0)] ^ <I>(r) d'^r, 


min 


(4.14) 


where the effective signal-to-noise ratio (snr) H is defined by 


H = A p2 $(0)/4>j^(0) = 
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+ 


-1 


(4.15) 


p2 WT (o 






1 + 


PB 
R ' 2 


4ttW 

A^BqO_ 


with = (AR)2/A the area of a resolution element in the object plane and 


o 



cp(u) d^u 


an area of the order of that of typical details in the object; Og^/Bo^ is the 
mean square contrast of the object, with Og^ = <P(0) . Table 1 gives the minimum 
relative mean square error for some simple object spectra. The more high 
frequencies attributed to the objects by <I>(r), the more slowly decreases 


with increasing snr H. 

The shot noise predominates over the effect of the field fluctuations 

when 


rg = A^BQ'a/47iW = A^Bq ' n/4trWhJ2 << 1, 


where n is the quantum efficiency of the photosensitive surface. If we con- 
sider a scene illuminated by moonlight on a clear night, we can put for Bq'/W 

12 

the value 1.6 x 10“3 watt*m~2y-l at A = 5150 A.U. This ratio rg then takes 
the value 0.8 x 10~^'' Pj which must be further diminished by the mean 
reflectivity of the scene. For illumination by full sunlight at the same 
wavelength the ratio rg is larger than this by a factor of about 10®. The 
quantum nature of light and its interaction with the recording medium is thus 
under most circumstances the principal hindrance to perfect image restoration. 
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5. Phase-Screen Turbulence 

When the object light passes through a turbulent medium, the point- 

spread function (psf) S 2 (r, u) appearing in Eq. (2.3) becomes a randomly time- 

varying function whose statistical properties derive from those of the 

turbulence. In order to assess the effect of the resulting distortion on the 

restorability of the received image, we take a simple model of the medium for 

which the average psf is readily calculated. 

The turbulence is supposed to be confined to a thin, homogeneous planar 

region between object and aperture, and its only effect is postulated as the 

introduction of random phase shifts into the rays passing through it. This 

is called the random phase ■. creen. As a result the spatial coherence function 

13 

of the object light at the erture acquires a factor 

p(r) = exp [- Y D(R'r/R) 3 , (5.1) 

where D(v) is the structure function of the random phase shifts 6(u) introduced 
at points u of the screen, 

D(v) = 2[0(O) - 0(v)3, 

0(v) = E[e(u) e(u + v)]. 

Here R*^ is the distance from object to phase screen. 

When in our image-processing system the integration time T is much 

longer than the time constant of the phase fluctuations, which in turn is 

much greater than the reciprocal bandwidth W"^ of the received and filtered 

( 2 ) 

light, the factor p(r) simply multiplies the optical transfer function p ^ (r) 
in the last integral of Eq. (4.9). If there are many independent phase 
screens lying between object and aperture and separated by many wavelengths 
of the light, the exponent in Eq. (5.1) contains a sum of structure functions 
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for the several screens; and if each screen introduces only a small, random 
phase shift, a good approximation to the net effect assigns to the factor 
p(r) a gaussian form, 

p(r) = exp (-r^/AL^), (5.2) 

where L is a characteristic length for the turbulence, whose distortions are 
the more severe, the smaller L is. Alternatively, this gaussian form can be 
taken as an approximation for a single phase screen. It corresponds to the 
light from each point of the object arriving with plane wavefronts tilted at 
random angles, which have a gaussian distribution about perpendicularity to 
the optical axis. 

The turbulence also affects the term in the mean square error, Eq. (4.9), 
containing the quadruple convolution I^^^(r, s) . This term arises from the 
product of four psf's 


Yi) S2*(r2> yP S2*(r3, V 2 ) S2(r^^, Y2~> 

occurring in the variance Var S(u) of the estimate. It must now be averaged 
over the ensemble of randomly fluctuating psf’s S 2 resulting from the turbu- 
lence. A tedious calculation making use of the inequalities W“^ << << T 


replaces I^^^(r, s) by a new function I^^^(r, s) involving the 


structure 


function of the random phase shifts. 


X 'ixp [r(|^ H - e)-<Ei - 52 - «>] + E> 

X d^rjd^r2d2t, (5.3) 

where R" is the distance from phase screen to aperture, R = R' + R", and 


,^(P. q) = exp {-[D(p) + D(q) - D(q - p) - D(q + p)]} . (5.4) 


20 



For a quadratic approximation to the structure function as in Eq. (5.2), 
D(p) “ p^, ^=1 and s) = I^^^(r, s) of Eq. (4.11). Furthermore, 

when the aperture is so large that diffraction imposes no significant limita- 
tions , 



s) = 




(5.5) 


The function S' can be shown to be at most equal to 1. Here we shall for 
simplicity set ^=1 and I^^^(r, s) = I^^^(r, s) , thus replacing the integral 
involving this kernel with its upper bound; the term involved is usually much 
smaller than the shot-noise term anyhow. 

We now assume that the diameter of the aperture is so much greater 
than the distortion length L that we can neglect diffraction of the image. 

The minimum mean square error is then, after suitable modification of Eq. 

(4.9) and subsequent equations. 


^ . = / [1 + Hlp(r)p $(r)/$(0)] ^ $(r) d^r, 

min ' _ ~ ~ ~ 


(5.6) 


where H is the snr defined in Eq. (4.15). As an example we assign to the 
spectral density of the radiance distribution the constant bandlimited form 


= Og^/ire^, |r| < e, 

= 0, |r| > e, (5.7) 

where with A the diameter of a typical detail in the object, e = XR/A. For 
the optical transfer function p(r) we use the gaussian form in Eq. (4.2). We 
then find for the minimum relative mean square error 


S ■ = 2(L/e) 2 £n { [H + exp (e^/2l2)]/(H + 1)}, 

min Jo 

which has been plotted versus H in Fig. 2 for various values of e/L = l/h, 
where i is the scale of the phase-screen turbulence projected onto the object 
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plane, % = XR/L. 

For comparison we present in Fig. 3 the minimum relative mean square 
error for the object spectral density 

$(r) = $(0) (r^ + (5.8) 

which extends to high spatial frequencies. The ratio ^ which was 

calculated numerically, decreases with H much more slowly than for the spectral 
density in Eq . (5.7). 
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Table 1 


Minimum Relative Mean Square Error 


Object Spectral Density 

S . 

min B 

$(r) = $(0), |rl < e 

(1 + H)-l 

= 0, jrl > e 


OCr) = 0(0) exp 

J-n (1 + H) 

0(r) = 0(0) (r^ + e^)"^ 

H ^ tan~^ (H^) 

0(r) = 0(0) (jr + e^)“^, n > 2 

(1 + 
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Figure Captions 


Fig. 1. 


Fig. 2., 


Fig. 3 


The object and the image-processing system: 0 - object plane, 

A = aperture plane, I = image plane. A narrowband temporal 
filter for object and background light is not shown. 

Minimum relative mean square erro. for restoring image dis- 
torted by a random phase screen. Object spectral density 
$(r) = 0(0), |rl < e; 0(r) = 0, jr] > e. Curves are indexed 
by the value of e/L = £/A. 

Minimum relative mean square error for restoring image dis- 
torted by a random phase screen. Object spectral density 
0(r) = 0(0) (r^ + e^)“^. Curves are indexed by the value of 

e/L = £/A. 
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